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Abstract 

A new set of nonlocal boundary conditions are proposed for the 
higher modes of the 3D inviscid primitive equations. Numerical schemes 
using the splitting-up method are proposed for these modes. Numeri- 
cal simulations of the full nonlinear primitive equations are performed 
on a nested set of domains, and the results are discussed. 

1 Introduction 

When the viscosity is present, the primitive equations have been the object 
of much attention, on the mathematical side. See the original articles [12, 
13], and the review articles about the mathematical theory of the PEs with 
viscosity appearing in [27] and in an updated form in [20]; see also the articles 
[1, 10, 11]. For the physical background on primitive equations, see e.g. [19] 
or [30]. In the absence of viscosity, little progress has been made on the 
analysis of the primitive equations since the negative result of Oliger and 
Sundstrom [18] showing that these equations are not well-posed for any set of 
local boundary conditions. However, the determination of suitable boundary 
conditions for the primitive equations is a very important problem for limited 
area models; see e.g. a discussion in [29]. 

In the broader context of limited-area numerical weather prediction mod- 
eling, the issue concerning the boundary conditions on the artificial bound- 
aries has been the focus of much research effort for decades. Since the 
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boundaries are artificial, the boundary conditions are expected to be of a 
transparent, or at least, non-reflecting type. It has been pointed out in the 
reference [18] cited above, and also in [6], that for many hyperbolic systems, 
truely non-reflecting boundary conditions will have to be non-local, i.e. in- 
volving states of the whole or part of the time period and/or the spatial 
domain. Generally speaking, non-local boundary conditions are difficult to 
implement in numerical simulations. Many authors have therefore proposed 
approximately non-reflecting boundary conditions. See, for example, [6], [7], 
[9], [16], and the references therein. These boundary conditions are also 
called absorbing boundary conditions, because they are designed to absorb 
the incident waves. Another approach that has been undertaken by some 
authors (see e.g. [17] and the references therein) is to introduce an absorbing 
layer, notably the perfectly matched layer (PML), surrounding the limited 
area. In this layer, the governing equations are modified to absorb any spuri- 
ous reflections. Our approach differs from those mentioned above in that we 
seek the truly non-reflecting boundary condtions, which are suitable for the 
governing equations in the sense of well-posedness, and are of a transparent 
type. As we have mentioned earlier, this type of boundary conditions will 
necessarily be non-local. They are valuable only if they can be shown to be 
practical for implementations in limited- area simulations. This is the main 
task of the current work. 

Following [26], two of the present authors (RT and JT) and A. Rou- 
sseau have investigated the inviscid primitive equations in space dimension 
two, and an infinite set of boundary conditions has been proposed. Well- 
posedness of the corresponding linearized equations has been established in 
[21] and numerical simulations have been performed in [22] for the linearized 
equations and for the full nonlinear equations. Note that the nonoccurrence 
of blow-up in the latter case supports the (yet unproved) conjecture that the 
proposed nonlocal boundary conditions are also suitable for the nonlinear 
PEs. 

Pursuing this approach, three of the present authors (QC, RT and JT), 
J. Laminie, and A. Rousseau considered a 2.5D model, with three orthogonal 
finite elements in the ^-direction, of the equations. The well-posedness result 
for the linearized equations was established in [2], and the numerical simu- 
lations of the nonlinear equations on a nested set of domains were discussed 
in [4]. 

The present article is related to the more theoretical ones [23] and [3] . In 
the first one [23], the authors obtained an infinite set of nonlocal boundary 
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conditions for the 3D inviscid primitive equations by studying the stationary 
problem associated with the linearized equations. In the second one [3], 
the authors gave due treatment of the special zero (barotropic) mode of the 
primitive equations and established the well-posedness of the corresponding 
linearized problem using the linear semi-group theory. Various numerical 
schemes through the projection method were also proposed, and the stability 
issue was studied for all of them. 

In the present work we intend to discuss the numerical simulations of 
the 3D nonlinear inviscid primitive equations on a nested set of domains. 
After performing the normal mode expansions of the unknowns, we are pre- 
sented with an infinite set of 2D equations. For the zero mode we use one 
of the schemes proposed in [3], which is semi-implicit, and is derived by the 
pressure-correction method. For the higher modes, i.e. the subcritical and 
supercritical modes, we use the splitting-up method for the discretizations 
and advance the unknowns along the x- and ^-directions in separate sub- 
steps. It then seems natural to impose boundary conditions by characteristics 
along the x- and ^/-directions separately. In the course of the article we recall 
the normal mode expansion leading to the infinite system of 2D equations 
(two spatial dimensions and time). Then we show how to discretize it in a 
form suitable for the implementation of the boundary conditions. 

Two simulations are performed. An initial simulation is carried out on a 
large domain with homogeneous boundary conditions. Using the data from 
the initial simulation as boundary conditions, we perform a second simulation 
of the same equations on a small interior domain. Then we compare these 
two results over the interior domain. 

Here the goals are twofold. On the one hand we want to numerically 
verify whether the boundary conditions, proven suitable for the linearized 
equations, are also suitable for the nonlinear equations. On the other hand, 
we want to numerically verify the transparency property of the proposed 
boundary conditions. Both goals are satisfactorily achieved. 

The article is organized as follows. In Section 2 we recall the 3D equations 
and their normal mode expansion. The issues of boundary conditions and 
well-posedness are also discussed. The numerical schemes are presented in 
Section 3. The settings for the numerical simulations, and the results of the 
simulations are discussed in Section 4. 
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2 The model 



The 3D primitive equations, linearized around a uniformly stratified flow (see 
[21], [22], [2] and [4]), read 



(2.1) 



v>t + U u x + (j) x - fv + B(u, v, w; u) = 0, 
vt + U v x + (j) y + fu + B(u, v, w\ v) + fU = 0, 
*l>t + Uo^ x + N 2 w + B(u, v, w] ij)) = 0, 
u x + v y + w z = 0, 



where u, v and w are the perturbation variables of the three velocity com- 
ponents, (j) is the perturbation variable of the pressure, ip is the perturbation 
variable of the temperature; / is the Coriolis force parameter, TV is the Brunt- 
Vaisala (buoyancy) frequency, assumed to be constant in the current study; 
B{u, v, w] 9) = u9 x + v9 y + w9 z for 9 = u,v, or ip. 

We will consider these equations in the domain M = M! x (— if, 0), 
M! = (0,Li) x (0,L 2 ), Li, L 2 , L 3 = H positive constants. Assuming flat 
bottom and the rigid lid hypothesis, we have 



(2.2) 



w 



at z = 0, -H. 



The boundary conditions for the other variables will be recalled and discussed 
below. 



2.1 Normal modes expansion 

Following [18] and [26], we consider the normal mode expansion of the solu- 
tions of the system (2.1). That is, we look for the solutions written in the 
following form: 



(2.3) 



(u,v,(f)) = ^2U n (z)(u n ,v n ,(f) n )(x,y,t), 

n>0 

(w, ip) = ^2 W n (z)(w n , ip n )(x, y, t). 



n>l 



Here U n and W n are solutions of the Sturm-Liouville boundary value problem 

^ = -\ 2 u(z), ze(-H,0), 
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respectively associated with the Neumann and Dirichlet boundary conditions. 
Therefore, we write the corresponding eigenfunctions as follows : 



/ Tin 

K a ■ 

yVn{z) = ^J^sm(X n z),U n (z) = cos ( X n z) , n > 1, 



(2.4) 



We then substitute the expressions (2.3) into (2.1), multiply each equation 
by U n (or W n for the 3rd and 5th equations), and integrate in z over the 
interval (— if, 0). We obtain the following systems: 



For n = 0, 
(2.5) 



du - du 

+ r 



dt 



dx dx 



~^~f Vo + J B(u,v,w;u)Uo(z)dz = 0, 



+ Uo^ + ^ + fu + J%(u,v,w;v)Uo(z)dz + fU ^ = 0, 

dx dy 
ipo = w = 0, 



For n > 1 
(2.6) 



at 

8v r . 



+ U 



fv n + J B(u,v,w;u)U n (z)dz = 0, 



+ t/o^r 1 + -7T- + fu n + / B(u,v,w;v)U n (z)dz = 0, 
, dt dx \ n dx dy 



+ / B(u, v, w- tP) W n (z) dz = 0. 

dy J-h 



The diagnostic unknowns (f> n and w n are given by 

1 



(2.7) 



A, 



and 

(2.8) w n = -^-{u nx + v ny ). 

2.2 Boundary conditions and well-posedness issues 

With the notation v = (u , ^o) T , the nonlinear equations of the zero mode 
can be written as 



(2.9) 
Here 

(2.10) Go 



(v t + U v x + fk x v + V0 O + G = 0, 
1 diw 0. 



f® H B(u, v, w] u) Uq(z) dz 
f® H B{u, v, w\ v) Uq(z) dz + fU Q y/H J ' 



and V and div are the 2D gradient and divergence operators, respectively. 
Without considering other modes, the nonlinear equations of the zero mode 
become 

v t + u v x + /kxnvni(^v)v = o, 

(2.11) ^ V H 

&\vv 0. 



where tp = (f) + fU VHy. 

We supplement the system (2.11) with the following boundary conditions: 



(2.12) 



uq = 0, at x = 0, Li, 

Vq = 0, at x = 0, and y — 0, L 2 . 



The well-posedness of the linearized system associated with (2.9), (2.12) has 
been studied in [3]. It is a standing conjecture that the boundary conditions 
(2.12) are also suitable for the nonlinear system (2.9), at least for a certain 
period of time. 

For the modes n > 1, we rewrite equation (2.6) in the matrix form as 
follows : 



Here 
(2.14) 

and 
(2.15) 




/ 


u 





-1 






A n 







Uo 







—N 2 






V 


A n 





Uo 



\ 



TP = 
5 ^ n 



/ o 








\ 






-1 












A n 






—N 2 






[° 


A n 





/ 



G n — 



-fv n + fl H B(u,v,w\u)U n (z)dz ^ 
fu n + fl H B(u,v,w]v)U n (z)dz 
f® H B{u, v, w\ ^)W n {z)dz ) 



We write 



(2.16) 



and 



(2.17) 




U r , 



N 




( 



v n + 



¥ 

Vn 

N } 



We define n c as the positive integer satisfying the following relations: 

n c 7r N (n c + 1)tt 
~W < U~ < H ' 

We will not study the non generic case where HN /nUo is an integer. We 
introduce the subcritical modes corresponding to 1 < n < n c , and the super- 
critical modes corresponding to n > n c . 

For the subcritical modes (n < n c ) we prescribe the following boundary con- 
ditions: 



(2.18) 



&(0,y,*)=0, 
v n (0,y,t) = 0, 
Vn(Li,y,t) = 0. 



(2.19) 



a n (x,L 2 ,t) = 0, 
/? n (x,0,t) = 0. 



And for the supercritical modes (n > n c ) we prescribe a slightly different set 
of boundary conditions: 



The well-posedness of the linearized system associated with (2.6) has been 
studied in [23] (see also [24]). The boundary conditions (2.18)-(2.21) are 
different from those proposed in [23] and [24]. We believe that the well- 
posedness of the linearized system corresponding to (2.6) and supplemented 
with the foregoing boundary conditions (2.18)-(2.21) can be established in 
the same way as in [23] and [24]; this problem will be studied elsewhere. We 
remark here that there are several sets of boundary conditions which make 
the linearized system well-posed. 

It is also a conjecture that the boundary conditions of [23] or [24] or the 
conditions (2.18)-(2.21) are suitable for the nonlinear equations for a certain 
time at least. 

3 The numerical schemes 

3.1 Numerical scheme for the zero mode 

Due to its resemblance with the classical Navier-Stokes equations and Euler 
equations, we discretize (2.9) by the pressure-correction method, which is a 
modified form of the classical projection method [15], [28], [8]. This modified 
form of the projection method is known to provide a better approximation of 
the pressure in the case of the Navier-Stokes equations and we choose to use it 
here, instead of the initial form of the projection method [5] , [25] . The bound- 
ary conditions (2.12) are different from those for either the Navier-Stokes 



(2.20) 




(2.21) 
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equations or the usual Euler equations; the pressure-correction method has 
to be adapted to the system (2.9). 

We let At = T/K, v k ~ v(x, y, kAt), and v k+ z represents an interme- 
diate value between v k and v k+1 , etc. At each step, the system is advanced 
in two substeps: 



(3.1) 

Here 
(3.2) 

and 
(3.3) 



At 

x=0 



G k 



+ U v h x +1 * +fkxv k + V(j) k + Gq -- 
= 0, 

f_ H B(u k , v k , w k ; u k ) Uq{z) dz 
f_ H B(u k , v k , w k ; v k ) Uq(z) dz + fU VH 



0, 



At 

V • v k+1 = 0, 
v k+1 • n = 0. 



+ V(^ +1 - 4) = 0, 



where n is the outer normal vector on dM! . 

It has been shown in [3] that if At, Ax and Ay satisfy the following 
conditions: 



(3.4) 



At 



< 



1 



(Ax 2 + Ay 2 ) 2 ~ c 2 K 4 ' 



At < 



1 



where c\ and K4 are constants independent of At, Ax and Ay, then the 
partially implicit scheme (3.1)-(3.3) is stable. For the details of the proof, 
and for the discussion of other related schemes, we refer the reader to [3]. 



3.2 Numerical scheme for the subcritical modes 

In this subsection we use the splitting method [14], [31], [25] to discretize 
(2.13), in the case of the subcritical modes. The supercritical case will be 
discussed in the next subsection. 

We have seen that the domain under consideration is (0, L\) x (0, L 2 ) x 
(0, T). We let /, J and K denote the numbers of grid pints in the x-direction, 
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the ^-direction and in time, and let Ax = Li/J, Ay = L 2 /J ) and At = 
T/K, denote the corresponding mesh. We denote the discrete grid points by 
(a;*, yj, tk) = (iAx, .7 Ay, fcAi). We let [/^ be the semi-discrete approximate 
value of {7 n at time and - the fully discrete approximate value of U n 
at (xi, i/j, tk). 

For each time step, two substeps are involved. The first substep is meant 
to advance the unknowns using only the advective terms in the x— direction 
and the zero order terms (the Coriolis force), that is, to solve the following 
semi-discrete equation: 

(35) U ^ ~ U " I E ^ +G -0 

(3.5) — +E n ^— + G n -0. 

The second substep is meant to advance the unknowns using only the ad- 
vective terms in the y— direction, that is, to solve the following semi-discrete 
equation: 

(3.6) L^ + ^.o. 

We now discuss the full discretization of the equations (3.5) and (3.6). 
We first apply the change of variables (2.16) to (3.5), and obtain 
(3-7) 

^ n + (Uo + y)y^~ = f<~ I B(u\v k ,w k ;u k )U n (z)dz 



At 








V n ~ 




At 








- 




At 



H 



+ N 



^ j\(u k ,v k ,w k ;^ k )W n (z)dz, 
+ Uo^~ = ~f^^ -J B(u\ v\ w k - v k )U n (z)dz, 



k+\ 



+ (U - yj^y = f<~ J B(u k , v k , w k - u k )U n {z)dz 

-y J° H B(u k ,v k ,w k ^ k )W n (z)dz, 
We recall that, for the subcritical modes, U — N/X n < 0. The up- wind 
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method applied to (3.7) yields the following fully discrete scheme: 
(3-8) 



' £^ + 2 _ tk 
^>n,i,j ^>n 

At " 

v . 1 — v ■ ■ 



7\ t 2 2 

™ \ ^n,i,j ~ Sn,i-l,j _ Q k,l • _ r , -i 



M + < e °+A:> a, 



fc+± 



.At 



Ax 



Vnjj ^n,i,j _|_ (jj^ _ N ^nj+ij Vn,i,j _ Q k,3 



At 



A n y Ax ^J'*- 1 ' ' J ' 

and j = 1, • • • , J + 1 in all cases, 



where 



/•0 

Sn'lj = + /<<j - / B(v* d , vfj, wl 3 -u\ 3 )U n (z)dz 

J-H 
1 f° 



(3.9) { 



S: 



k,2 



v k ■ ■ - f 



-f 

J-H 



B{ul J) vl J) wl J ' ) vl J )U n {z)dz ) 



Snlj = Vn,ij + /<<J " / B{u k ^ V k hJ) W^ J )U n (z)dz 

J-H 

1 f° 



Jc+± k+i 



k+i 



The boundary conditions for t; n 2 , v n 2 and rj n 2 are, for < j < J, 
(3.10) 

We then apply the change of variables (2.17) to (3.6) and obtain 



e^H~ 2 n 2 n kJ[ ~ 2 n 

n,0,j — U > V n,0,j — U ' Vnjj ~ U> 



(3.11) 



- Un +2 



At 



0, 



a k + l - an 2 N da k + l 



At 



A n dy 



0, 



At 



+ 



N = Q 
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Applying the up-wind method to the system (3.11) yields 



(3.12) 



11 • • II 



At 



n,i,j n,i,j iV n,i,j 



Jc+1 



At 



A. 



Ay 



= 0, 



,3/C+l _ O- ' 2 AT /Ofc+1 _ OK+1 
^n,i,j ' n,i,j _|_ ^ n i*J ' n ,i^3~ 1 q 



At A n Ay 

The boundary conditions for ce* +1 , are 



(3.13) 



/c+1 



o, 



for < j < J, 



^ = 0, for < i < I. 



We remark here that u^ +1 does not need any boundary conditions. 



3.3 Numerical scheme for the supercritical modes 

The fully discrete numerical schemes for the supercritical modes can be de- 
rived by the same approach presented in the previous subsection. The results 
for the supercritical modes are similar to those for the subcritical modes, and 
are simpler because all the eigenvalues of the coefficient matrix E n are pos- 
itive. We shall omit the intermediate details, and present the numerical 
schemes directly. Only the differences with those for the subcritical modes 
will be pointed out. 

As for the subcritical modes, the numerical schemes for the supercritical 
modes also involve two substeps. The first substep consists of the following 
scheme: 
(3-14) 

~~Af + {Uo + T n ) —Ai -S n , iJt i-2,.~ ,1 + 1, 

At + ° Ar " 1 ~ 2 ' ' ' ' ' 1 + 



^H - 2 a T k~\- 2 ^"^2 

^n,i,j ^n,ij /f T ^ \ Vn,i,j Vn,i—l,j nk,3 

At + " YJ a^ - S ™ 

and % = 1, • • • , J + 1 in all cases , 
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Here > S nlji Snh> aild S nh ^ defined in 3 ' 9 ' We n ° te here that dp^/dx 

is discretized differently in (3.8) and in (3.14), due to the fact that Uq — 1/X n 
has different signs in the sub- and super-critical modes. The boundary 



k+\ 



conditions for ^ n 2 , v n 2 and rj n 2 are, for < j < J, 



(3.15) 



The second substep consists of the following scheme: 



(3.16) 



At 

n,i,j n,x,j iV ^n^j+l n,i,j 
"n,i,j "n,i,j _|_ " n ,h3 q 



= o, 



At X n Ay 

The boundary conditions for a^ +1 , /3^ +1 are 



(3.17) 



OL 



/c+1 
Dk+1 



for < j < J, 



TO = °> for < z < /. 



3.4 Treatment of the integral of the nonlinear term 

In this section, we will deal with the integral of the nonlinear term. There 
are five kinds of integrals to be considered. We have the following lemma. 

Lemma 3.1 Assume that u,v,(j),w, and ip have the expressions (2.3) and 
U n and W n for n > are as defined in (2.4). Then 



(3.18) 

r0 



B{u, v, w] u)Uo(z)dz 



H 



=—y 



du m du m 



m>0 



dx 



dy 



m>l 



(3.19) 

r0 



H 



m>0 



dx 



dy 



m>l 
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(3.20) 



/o 1 n 

B(u,v,w;u)U n (z)dz = —== V( 

_ v^/ oum du m 1 y- 



'2H 



du m 
dx 

du m 



du m 
dy 

du m 



m=n 
oo 



2// ^ dx dy 

m>0 



) 



rn W rn J rn U m 



m>n+l 



m>l 



1 n 



m>l 



(3.21) 



/o 1 " 

OO 



n 



<9u 



m 



m—n 



dx 



dv, 



m 



m—n 



v 1 9?; m 5^ m . 



m=n 

oo 



ro>0 



= E * 



'2# 



ra>n+l 



^ oo 1 72 



m>l 



m>l 



(3.22) 



f B(u,v,w;xp)W n (z)dz = — = V7?x m _ n ^^ + v m _J^^) 
J-h V2H ^ dx dy 



m—n 

oo 



1 v^/ d ^m , d*l) m 1 v^/ d ^rn , ^ 

-= ) (u n _ m — \-v n - m —— ) -= ) [u, 

/ 2H ^ dx dy x/2H ^ 

m=l v 



n-\-m ' 



m=l 



dx 



+ v n+r, 



dy 



1 n 

OH ^ 



f 2H 



^ OO ^ oo 

2H . , v2xi ^ 



m>l m>l m>n 

Lemma 3.1 can be verified by direct calculations. 

Remark 3.2 in large-scale GFD simulations, in which a large number of 
modes are involved, the preceding convolution products would be too costly in 
terms of CPU time to be appropriate. To avoid them, it is then necessary to 
transform the Fourier coefficients u n) etc., back into the physical space, com- 
pute the nonlinear products in the physical space, and calculate the integrals 
on the left side of (3. 18) -(3. 22). In our study, only a small number (< 10) 
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of modes are considered, and thus the formulas (3.18)-(3.22) are appropriate 
and sufficient. 

4 Numerical simulations in a nested environ- 
ment 

Two different simulations are performed. The first one is carried out on 
the larger domain Ai = (0,Li) x (0, L 2 ) x (—77,0) (see Figure 1), and a 
set of homogeneous boundary conditions prescribed at (x,y) G d M! where 
M! = (0, Li) x (0, L 2 ). The simulations will be described and the results will 
be presented in details in Section 4.1. The data obtained through this simu- 
lation will provide the nonhomogeneous boundary conditions for the second 
simulation on the middle half domain, denoted by A4i = (Li/4, 3Li/4) x 
(L 2 /4,3L 2 /4) x (-if, 0) (see also Figure 1), of M. This simulation will be 
described and the numerical results will be presented in detail in Section 4.2. 

In Section 4.3, the numerical results from these two simulations are then 
compared, and the coincidence of the numerical results demonstrates the 
transparent properties of the proposed boundary conditions, and supports 
the conjecture of their suitability for the nonlinear equations. 

The physical parameters that we used in the simulations are the following 
ones: L\ = 1000km, L 2 = 500km, H = 10km. We take the constant reference 
velocity f/ = 20 m/s, the Coriolis parameter / = 10 -4 , and the Brunt- 
Vaisala (buoyancy) frequency TV = 10 -2 . The final time for the simulations 
is T = 5 x 10 4 s, and we take 1600 time steps. In the vertical direction we 
take 40 segments. In the computations, we will deal with iV max = 5 (the 
number of modes), which is sufficient from the physical point of view. 
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fi larger domain with homogeneous BC's 



fii middle half domain with non-homog BC's 

Figure 1: The larger domain Ai and the middle half domain Ai\ 
4.1 Simulation on the larger domain 

In the simulation, the initial conditions are given for these scalar functions: 
(4.1) 

x 2ir . 2nx 2ny . Anx Airy nz 
u(x,y,z,0) = — — sm( — )cos( — ) + sm( — )cos( — )cos( — ), 

1j\ Lio L\ Li2 Li\ Li2 n. 

— if 27TX 27TX 27TX \ ,2iry 
v(x, y, z, 0) = — I sin ( — ) + — cos ( — ) I sin (— ) 

Z/2 / o Anx Anx Any ,^ z \ 

+ T 1 \ + sm (^) sm (^) cos (#) 

/ ^\ —4:H . Airx , A^ x \\ A' k V\ ■ i vz \ 

w(x, y, z, 0) = — — (sin (— - ) + cos (-— )) cos (— ) sin (— ), 

±j\ ±j\ ±j\ L/2 £1 

- 2ttx ,27n/ w ttz /2ttz 
</){x,y,z,0) = C/ sin( — )sin(— )(cos(— ) - cos(— )), 

^(x,?/,z,0) = — sm (— ) sm ( — )(2 sm (— ) - sm (-)). 

We note here that these initial functions u, v, w, 0, and i/j satisfy the homoge- 
neous boundary conditions for each mode n > 0. Specifically, for the zeroth 
mode, i.e. when n = 0, 



(4.2) 



«o(0,y,*) = 0, uo{L u y,t) = 0, 

•u o (0,y,£) = 0, w (x,0,t) = 0, v (x,L 2 ,t) = 0; 
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Figure 2: The initial state of velocity field in the larger domain Ai 



f n (0,y,t) = 0, v n (0,y,t) = 0, r] n (L u y,t) = 0, 
a n (x, L 2 ,t) = 0, p n (x,0,t) = 0; 



for the subcritical modes, i.e. when 1 < n < n 
(4.3) 

and for the supercritical modes, i.e. when n> n C) 
(4.4) 



f n (0, y, t) = 0, v n (0, y, t) = 0, 7? n (0, y, t) = 0, 
a n {x, L 2 , t) = 0, /3 n (x, 0, i) = 0. 



In this simulation, we take 400 segments in the x-direction, and 200 seg- 
ments in the y direction. When restricted to the middle half domain, the 
functions (4.1) also provide the initial conditions for the simulations on the 
middle half domain. 

The simulation results over the larger domain A4 are plotted in Figures 
2 to 17. Figure 2 is the cone plot with isosurface of the initial state of the 
velocity field, and Figures 3 and 4 are the slice-plane plots of the initial 
state of (j) and ip. Figures 5 to 9 are the contour plots of u ) v, w, if) and 0, 
respectively, on the plane z = — 2, 500m, at t — 0. 

Figure 10 is the cone plot with isosurface of the velocity field at the final 
time t — T \ and Figures 11 and 12 are the slice-plane plots of the state of 
and i/j at the final time t = T. Figures 13 to 17 are the contour plots of u, v, 
w, i/j and 0, respectively, on the plane z = —2, 500m, at t — T. 
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Figure 3: The initial state of ip in the larger domain A4 
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Figure 4: The initial state of (j) in the larger domain A4 
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Figure 11: The state of ip in the larger domain Ai at t = T. 
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Figure 12: The state of (j) in the larger domain A4 at t = T. 
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Figure 13: Contour plot of u at z = —2500m, at t — T. 
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Figure 17: Contour plot of (j) at z = —2500m, at t — T. 
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4.2 Simulation in the middle-half domain 



As explained in the Introduction, we next do simulations on the middle half 
domain Mr of M, Mr = M[ x (-#, 0), M\ = ^) x (^f, The 
boundary values of the unknown functions u, v, w, 0, andip are inferred from 
the previous simulation. More specifically, the boundary conditions are, for 
the zeroth mode (n = 0), 



(4.5) 



«o(W 4 > Vj,tk) = «o(^i/4, *fc), 
uo(3.Li/4, = ii{,(3Li/4, 

voiLx/^yj, t k ) = ^(Li/4,%-,4), 
vofo, L 2 /4, = ^(^, L 2 /4, 
^o(^,3L 2 /4,t fc ) = ^(xi,3L 2 /4,* fc ). 



For the sub critical modes (1 < n < n c ), 



(4.6) 



£„(Li/4, f fc ) = ^(Li/4, * fc ), 
u„(Li/4, t fe ) = ^(Li/4, 4), 
?7 n (3Li/4, ^-,4) = ^(3Li/4, 
a ra (:r;,3L 2 /4,4) = a4(^, 3L 2 /4, t k ), 
P n (x it L 2 /4, t fc ) = p l n ( Xi , L 2 /4, 4), 



and for the supercritical modes (n > n c 



(4.7) 



£„(Li/4, f fc ) = &(Li/4, t fc ), 
u„(Li/4, t fc ) = u^(Li/4, t fc ), 
Vn(Li/4,yj,t k ) = 77^(^/4,^,4), 
a ra (a;i,3L 2 /4,4) = a;^, 3L 2 /4, £ fc ), 
/3 n (^, L 2 /4, * fc ) = L 2 /4, * fc ). 



In the above, Xj, and tk denote the discrete grid points in space and time. 
The superscript / denotes the previous simulation in the larger domain M. 
In this simulation, we take 200 segments in the x-direction, and 100 segments 
in the y direction. 

The simulation results over the middle half domain Mi are plotted in 
Figures 18 to 25. Figure 18 is the cone plot with isosurface of the velocity 
field in the middle half domain Mi at the final time t — T, and Figures 19 
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Figure 18: The velocity field with cone plot in the middle half domain Ai\ 
at t = T 



and 20 are the slice-plane plots of the state of and ip in the middle half 
domain A4± at the final time t = T. Figures 21 and 25 are the contour plots 
of u, v, ijj and 0, respectively, on the plane z = —2, 500m restricted to the 
middle half domain .Mi, at t = T. 



4.3 Comparison 

In this subsection, we compare these two distinct simulations namely the 
results of the simulations on the larger domain A4 restricted to the middle 
half domain A4i and the results of the simulations on the middle half domain 
Aii, as obtained by the second simulation above. 

Let u ext , v ext , w ext , (f) ext and vp ext be the numerical approximations of 
the variables u, v, 0, and %j) on the larger domain Ai, respectively, and 
u int , v int , w int , (j) int and ip znt be the numerical approximations of the variables 
u, v, w, and ijj on the middle half domain A^i, respectively. In Figures 
26-30, we plot the evolution of the unknowns, and of their relative errors (see 
below), in both the L 2 and L°° norms. The relative errors are defined as 



\U int — U ext \M 1 \\LP 
\\u ext \\LP 



etc. where p = 2, oc. 



We observe that both the L 2 and the L°° norms of the prognostic variables 
u ) v and tj) are diminishing in time. This can be explained by the homoge- 
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X 

Figure 19: The state of ip in the middle half domain A4± at t — T. 
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Figure 20: The state of 4> in the middle half domain A4± at t = T. 
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Figure 21: Contour plot of u at z = —2500m, at t — T. 

neous boundary conditions imposed on the boundary of the larger domain 
Ai and by the fact that the velocity field is flowing out of the domain to 
the right, with a constant velocity U$. The relative errors of the prognostic 
variables, in both the L 2 and the L°° norms, are of the magnitude O(10 -2 ) 
or smaller, which means that results on the larger domain A4 and the results 
on the middle half domain A4[ match very well. 

For the diagnostic variables, the L 2 and the L°° norms of w are also 
diminishing in time, as for the prognostic variables; the relative errors for w 
are large (O(10 -1 )) as compared to those for the prognostic variables, but 
they are still well controlled. The bizarre behavior of the L 2 and the L°° 
norms of (j) can be explained by the absence of an evolution equation for 
(j) and the lack of natural boundary conditions for O in (3.3). The relative 
errors for (j) are however very well controlled. Note that the relative errors 
for v, and ^, in both the L 2 and L°° norms, are of the order of 
O(10 2 ), and the relative errors for w are of the order of O(10 _1 ). 

A graph of the absolute divergence averaged over the guest integration 



36 




Figure 22: Contour plot of v at z = —2500m, at t — T. 



area is presented in Figure 31 for three cases: larger domain, the middle-half 
domain from direct computation, and the middle-half domain using the data 
from the larger domain. We observe that the mean absolute divergence for 
three cases are small and diminishing in time. This can be explained by 
the divergence free condition employed on the proposed numerical schemes. 
Furthermore, Figure 31 shows that the behaviors of the mean absolute di- 
vergence match well on the middle-half domain J\d\. 



5 Conclusions 

In conclusion, the absence of blowing up demonstrates that the boundary 
conditions proposed in Section 2.2 are suitable for the problem, and the 
numerical scheme proposed in Section 3 is stable. The fact that the numerical 
results match very well on the middle half domain A4\ demonstrates the 
transparency property of the boundary conditions. 
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Figure 23: Contour plot of w at z = —2500m, at t — T. 

From the results of this idealized model it is straightforward to outline 
the algorithmic path to be taken in the application of this method to the 
full primitive equations. The simplest approach would be first to re- write 
the model equations so that they are formally equivalent to the system (2.1). 
This involves specifying a reasonable, local mean stratification TV 2 , and mean 
zonal wind, Uq. 

Next a vertical mode decomposition is performed to identify the sub- 
critical/supercritical mode division. Next, the appropriate lateral boundary 
conditions are applied. Lastly the modal decomposition is summed to re- 
construct the boundary values of the field variables. As necessary, the local 
mean stratification and zonal wind can be adjusted. 

This would be superior to methods that absorb wave energy through 
nudging since the artificial damping also inevitably causes the interior solu- 
tion to decay and the sponge layer, itself, induces wave reflection. 
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Figure 24: Contour plot of ip at z = —2500m, at t — T. 
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